function [Tn, cv] = as10_test_corr(m, g_ar, N_mkt, rand_seed, reps, sig_level, ...
    correlated_fc, N_dev2, corr_m, corr_g_ar)
% AS test statistics and critical values

seed = RandStream('mt19937ar','Seed',rand_seed);

small_cons = 1e-10;
r0 = length(g_ar);

[K, nm] = size(m);% N_mkt*N and # moment functions (if FY) and N_mkt  and # moment functions (if CT or AT)
if K ~= N_mkt; N = K/N_mkt; end

count_r = 0; % number of g functions
for r = 1 : r0
    count_r = count_r + size(g_ar{r}, 2);
end

% interact the moments with instruments
M = nan(K, nm, count_r);
count_rm = 0;
for r = 1 : r0
    for a = 1 : size(g_ar{r}, 2)
        count_rm = count_rm + 1;
        M(:,:,count_rm) = m.*g_ar{r}(:, a);
    end
end

if correlated_fc
    n_corr_iv = size(corr_g_ar, 2);
    corr_miv = nan(N_mkt, 2, n_corr_iv);

    for a = 1 : n_corr_iv
        % interact the moments with instruments
        corr_miv_a = corr_m.*corr_g_ar(:, a);

        % average across potential products within a market
        corr_miv(:, :, a) = mean(permute(reshape(corr_miv_a', [2, N_dev2 N_mkt]), [3 1 2]), 3); % N_mkt x 2
    end
else
    corr_miv = [];
end
clear m g_ar

if K ~= N_mkt
    if correlated_fc
        M = [squeeze(mean(reshape(M, [N_mkt, N, nm*count_r]), 2)),...
            reshape(corr_miv, [N_mkt, 2*n_corr_iv])];% N_mkt x (nm*count_r + 2*n_corr_iv)
    else
        M = squeeze(mean(reshape(M, [N_mkt, N, nm*count_r]), 2));% N_mkt x (nm*count_r)
    end
else
    M = reshape(M, [N_mkt, nm*count_r]);
end

N_moments = size(M, 2);

%% test statistic
mean_M = mean(M); % 1 x N_moments

cov_M = cov(M);

Tn = sum((sqrt(N_mkt)*mean_M.*(mean_M>0)./sqrt(diag(cov_M) + small_cons)').^2);

%% Critical Value
Dinvsqrt = 1./sqrt(diag(cov_M)+small_cons);

Omg = diag(Dinvsqrt)*cov_M*diag(Dinvsqrt);
Omg = Omg*1/2 + Omg'*1/2 + eye(N_moments)*(small_cons);

% GMS selection
xin = 1/sqrt(log(N_mkt))*Dinvsqrt'.*mean_M.*(mean_M<0).*sqrt(N_mkt); % 1 x N_moments

sim_Z = randn(seed, reps, N_moments)*chol(Omg) + xin; % reps x N_moments

Tn_S2_reps = sum(((sim_Z>0).*sim_Z./sqrt(diag(Omg)')).^2, 2);

cv = quantile(Tn_S2_reps, 1 - sig_level);








